// 针对QHZX-FY3D-NDVI2 产品开发的 NDVI 数据拼接合成程序，
// 支持轨道数据合成到日产品，支持日产品合成多天产品
// 合成算法需要观测角度
// 合成逻辑 查找VZ最小的两个值，然后保留最大的一个NDVI，同时保留该VZ
// 既然是NDVI那么只做第一个波段
// 2022-3-14
// v1.2.0 使用最大值合成，不考虑vza

#include <fstream>
#include <sstream>
#include <vector>
#include <string>
#include "wstringutils.h"
#include "wGdalRaster.h"
#define WGDALRASTER_H_IMPLEMENTATION
#include "wGdalRaster.h"
#include "wdirtools.h"
#include <iostream>
#include "gdal_priv.h"
#include "ogr_spatialref.h"
#include <memory>


using namespace std;

struct LLValue{
	double x,y,ndvi,vza ;
};

struct TwoFile{
	TwoFile(){}
	TwoFile(string nfilename,string vfilename):ndvifile(nfilename),vzafile(vfilename){}
	string ndvifile,vzafile ;
} ;

bool ReadFileList(string listfilename ,vector<TwoFile>& filenameVec)
{
	std::ifstream infs(listfilename.c_str());
	if( infs.good()==false ) return false ;
	std::string line;
	while (std::getline(infs, line))
	{
		if( line.length() > 2 )
		{
			TwoFile twofile ;
			vector<string> tstrvec = wStringUtils::splitString(line,",") ; 
			if(tstrvec.size() !=2 ){
				cout<<"ReadFileList exception: tstrvec not 2."<<endl ;
				cout<<"line:"<< line <<endl ;
				cout<<"line.len:"<< line.length() <<endl ;
				return false ;
			}
			twofile.ndvifile = tstrvec[0] ;
			twofile.vzafile = tstrvec[1] ;
			filenameVec.push_back( twofile ) ;
		}
	}
	return true ;
}


bool FillOneFileIntoOutput( 
	wGdalRaster* outRas0 ,
	string ndvifile ,
	double ndvifilldata,
	int iinfile_for_progress,
	int numinfile_for_progress,
	string& error
	)
{
	shared_ptr<wGdalRaster> inNdviRas(wGdalRasterFactory::OpenFile( ndvifile )) ;
	if( inNdviRas.get() == 0 ){
		error = ndvifile + " is bad";
		return false ;
	}
	const double* inTrans  = inNdviRas->getTrans() ;//maybe not wgs84
	char* inProj = inNdviRas->getProj() ;
	const double* outTrans = outRas0->getTrans() ;
	char* outProj = outRas0->getProj() ;
	OGRSpatialReference outSrs(outProj);
	OGRSpatialReference inSrs(inProj);

	//OGRCoordinateTransformation* pTransformer=OGRCreateCoordinateTransformation(&inSrs,
	//	&outSrs);
	OGRCoordinateTransformation* pTransformerO2I =
		OGRCreateCoordinateTransformation(
			&outSrs,
			&inSrs);

	const int inxsize = inNdviRas->getXSize() ;
	const int inysize = inNdviRas->getYSize() ;
	const int outxsize = outRas0->getXSize() ;
	const int outysize = outRas0->getYSize() ;

	int fileper0 = iinfile_for_progress*100.0 / numinfile_for_progress ;
	int per0 =  -1;
	for(int oy = 0 ; oy < outysize ; ++ oy )
	{
		double olatcenter = outTrans[3] + outTrans[5] * oy + outTrans[5]/2;
		for(int ox =0 ; ox < outxsize ; ++ ox )
		{
			double oloncenter = outTrans[0] + outTrans[1] * ox + outTrans[1]/2 ;
			double incx = oloncenter ;
			double incy = olatcenter ;
			pTransformerO2I->Transform(1, &incx, &incy ) ;
			incx = incx - inTrans[1]/2 ;
			incy = incy - inTrans[5]/2 ;

			int iix = (incx - inTrans[0])/inTrans[1] ;
			int iiy = (incy - inTrans[3])/inTrans[5] ;

			if( iix >=0 && iix < inxsize && iiy >=0 && iiy < inysize )
			{
				double newval = inNdviRas->getValued(iix,iiy,0) ;
				if( newval != ndvifilldata )
				{
					if( outRas0->getValued(ox,oy,0) == ndvifilldata ){
						outRas0->setValued(ox,oy,0,newval) ;
					}
					else if( newval > outRas0->getValued(ox,oy,0) ) {
						outRas0->setValued(ox,oy,0,newval) ;
					}
				}
			}			
		}
		int per1 = fileper0 + oy * 100.0 / outysize / numinfile_for_progress ;
		if( per1 != per0 ){
			per0 = per1 ;
			cout<<per0<<endl ;
		}
	}

	OCTDestroyCoordinateTransformation(pTransformerO2I);
	return true ;
}


int main(int argc , char** argv )
{
	cout<<"A program to composite multi NDVI product with VZA data, using VCMVC method."<<endl ;
	cout<<  "v1.0 2022-3-14"<<endl ;
	cout<<"v1.0.1 2022-3-16"<<endl ;
	cout<<"v1.1.0 2022-3-17 use deflate compress."<<endl ;
	cout<<"v1.2.0 2022-3-22 from output to transform, use MVC without VZA."<<endl ;
	cout<<"version: v1.3.0.1 2024-3-3 use dir prefix tail to scan input files."<<endl;
	cout<<"usage-1:"<<endl;
	cout<<"  qh5_ndvi_vcmvc infilelist.txt inoutfilldata invzafilldata out_xmin out_ymin out_xmax out_ymax out_reso output_ndvifile.tif output_vzafile.tif"<<endl ;
	cout<<"infilelist.txt: one line has: ndvi_filename,vza_filename, use , to separate."<<endl ;
	cout<<"usage-2:"<<endl;
	cout<<"  qh5_ndvi_vcmvc mode /indir/ prefix tail filldata outxmin outymin outxmax outymax outreso outdatafile"<<endl;
	cout<<"mode should use 2."<<endl;
	cout<<"e.g.:"<<endl ;
	cout<<"/some/ndvidir/ndvi1.tif,/some/vzadir/vza1.tif"<<endl ;
	cout<<"/some/ndvidir/ndvi2.tif,/some/vzadir/vza2.tif"<<endl ;
	cout<<"/some/ndvidir/ndvi3.tif,/some/vzadir/vza3.tif"<<endl ;

	GDALAllRegister() ;
	if( argc==11 || argc==12 ) {
		//11 for useage-1
		//12 for useage-2
		//ok
	}else{
		cout<<"argc neither 11 nor 12. argc="<<argc <<endl ;
		return 11 ;
	}

	vector<TwoFile> infileVec ;
	double ndvifilldata = 0;
	double vzaafilldata = 0;//deprecated.
	double outxmin = 0 ;
	double outymin = 0 ;
	double outxmax = 0 ;
	double outymax = 0 ;
	double outreso = 0 ;
	string outndvifile  ;
	string outvzaafile  ;//deprecated

	if( argc==11 )
	{
		string inlistfile = argv[1] ;
		ndvifilldata = atof(argv[2]) ;
		vzaafilldata = atof(argv[3]) ;
		outxmin = atof(argv[4]) ;
		outymin = atof(argv[5]) ;
		outxmax = atof(argv[6]) ;
		outymax = atof(argv[7]) ;
		outreso = atof(argv[8]) ;
		outndvifile = argv[9] ;
		outvzaafile = argv[10] ;

		cout<<"------------"<<endl ;
		cout<<"inlistfile :"<<inlistfile<<endl ;
		cout<<"ndvifilldata :"<<ndvifilldata<<endl ;
		cout<<"vzaafilldata :"<<vzaafilldata<<endl ;
		cout<<"outxmin :"<<outxmin<<endl ;
		cout<<"outymin :"<<outymin<<endl ;
		cout<<"outxmax :"<<outxmax<<endl ;
		cout<<"outymax :"<<outymax<<endl ;
		cout<<"outreso :"<<outreso<<endl ;
		cout<<"outndvifile :"<<outndvifile<<endl ;
		cout<<"outvzaafile :"<<outvzaafile<<endl ;//deprecated
		cout<<"------------"<<endl ;
		bool listok = ReadFileList( inlistfile , infileVec ) ;
		if( listok==false ){
			cout<<"failed to read inlistfile:"<<inlistfile<<endl ;
			return 12 ;
		}
		if( infileVec.size()==0 ){
			cout<<"infileVec size equal 0. "<<endl ;
			return 13 ;
		}
	}else if( argc==12 )
	{//for version v1.3.0
		string mode = argv[1] ;//always be 2, useless
		string indir = argv[2] ;
		string prefix = argv[3];
		string tail = argv[4] ;
		ndvifilldata = atof(argv[5]) ;
		outxmin = atof(argv[6]) ;
		outymin = atof(argv[7]) ;
		outxmax = atof(argv[8]) ;
		outymax = atof(argv[9]) ;
		outreso = atof(argv[10]) ;
		outndvifile = argv[11] ;

		cout<<"------------"<<endl ;
		cout<<"mode :"<<mode<<endl ;
		cout<<"indir :"<<indir<<endl ;
		cout<<"prefix :"<<prefix<<endl ;
		cout<<"tail :"<<tail<<endl ;
		cout<<"ndvifilldata :"<<ndvifilldata<<endl ;
		cout<<"outxmin :"<<outxmin<<endl ;
		cout<<"outymin :"<<outymin<<endl ;
		cout<<"outxmax :"<<outxmax<<endl ;
		cout<<"outymax :"<<outymax<<endl ;
		cout<<"outreso :"<<outreso<<endl ;
		cout<<"outndvifile :"<<outndvifile<<endl ;
		cout<<"------------"<<endl ;
		vector<string> vec1 ;
		wDirTools::getDirFilesRecursive2(
			indir, prefix, tail, vec1) ;
		cout<<"scan files count : "<<vec1.size()<<endl;
		if( vec1.size()==0 ){
			cout<<"vec1 size equal 0. "<<endl ;
			return 13 ;
		}
		for(int it =0;it<vec1.size();++it ) infileVec.push_back(TwoFile(vec1[it],"")) ;
	}

	const int outxsize = (outxmax-outxmin)/outreso + 1 ;
	const int outysize = (outymax-outymin)/outreso + 1 ;
	cout<<"outsize "<<outxsize<<","<<outysize<<endl;
	if( outxsize<=1 || outysize<=1 ){
		cout<<"output size invalid "<<outxsize<<","<<outysize<<endl;
		return 14 ;
	}

	GDALDataType thedatatype = GDT_Byte ;
	// only one band ok.int theBandNum = 0 ;
	{
		shared_ptr<wGdalRaster> tempPtr(wGdalRasterFactory::OpenFile( infileVec[0].ndvifile )) ;
		if( tempPtr.get() == 0 ){
			cout<<"bad ndvifile 0."<<endl ;
			return 14 ;
		}
		thedatatype = tempPtr->getDataType() ;
	}

	{
		shared_ptr<wGdalRaster> outRas0( wGdalRasterFactory::Create(outxsize,outysize,1,thedatatype) ) ;
		if( outRas0.get() == 0 )
		{
			cout<<"create outRas0 failed."<<endl ;
			return false ;
		}
		outRas0->fill(0 , ndvifilldata);
		outRas0->copyProj("GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]"); 
		const double outTrans[6] = {outxmin , outreso , 0.0 , outymax , 0.0 , -outreso } ;
		outRas0->copyTrans(outTrans) ;

		for(int iin = 0 ; iin < infileVec.size() ; ++ iin )
		{
			cout<<"ifile "<<iin<<"/"<<infileVec.size()<<endl ;
			string error1;
			bool ok1 = FillOneFileIntoOutput( 
				outRas0.get(),
				infileVec[iin].ndvifile,
				ndvifilldata,
				iin,
				infileVec.size(),
				error1
				);
			if( ok1==false ){
				cout<<"bad fill for "<<infileVec[iin].ndvifile<<","<<infileVec[iin].vzafile<<endl ;
				cout<<"error:"<<error1<<endl ;
				return 16 ;
			}
		}
		cout<<"100."<<endl ;
		cout<<"saving..."<<endl ;
		bool saveok1 = outRas0->save(outndvifile,ndvifilldata,true) ;
		if( saveok1==false ){
			cout<<"failed to save "<<outndvifile<<endl ;
			return 17 ;
		}
	}
	cout<<"success"<<endl ;

	return 0 ;
}